Test Problem

The descritization of the PDE on one by one squre with n interior ponts in z direction and m interior points in v direction.

{n, m}  =  {20, 10};
{hz, hv}   =  {(1/(n + 1)), (1/(m + 1))} // 
          N;
{aa, bb, c, d, 
          e}  =  {-((2/(hz^2)) + (2/(hv^2)) - 10), (1/(hz^2)) + 10 hz, 
            (1/(hz^2)) - 10 hz, (1/(hv^2)), (1/(hv^2))} // N;

A  =  
        SparseArray[{{i_, i_} -> 
              aa, {i_, 
                  j_} /; (i - j == 
                      1) && (! IntegerQ[j/m]) -> 
              bb, {i_, 
                  j_} /; (i - j == -1) && (! IntegerQ[
                        i/m]) -> 
              c, {i_, j_} /; (i - j == m) -> 
              d, {i_, j_} /; (i - j == -m) -> e}, {n*m, 
            n*m}];
B  =  Table[Random[], {n*m}, {1}];
CC  =  Transpose[B];

Structure of A

ArrayPlot[Take[A, 50, 50],  AspectRatio -> Automatic]

Solution of Lyapunov equation

A X  + X (A^T) =  -B (B^T) 
(A^T) X  + X A =  -(C^T) C

<< LargeScaleControl`LyapunovSparse`

Xc  =  LowRankLyapunovSolve[A, B,  Tolerance -> (10^(-10))]; // 
      AbsoluteTiming

nrm  =  
        Norm[A . Xc + Xc . Transpose[A] + 
            B . Transpose[B]];
If[Head[nrm]  ==  Real,  nrm]

Xo  =  LowRankLyapunovSolve[A, CC, Tolerance -> (10^(-10))]; // 
      AbsoluteTiming

nrm  =  
        Norm[Transpose[A] . Xo  +  Xo . A + 
            Transpose[CC] . CC];
If[Head[nrm]  ==  Real, nrm]

Xfull  =  LyapunovSolve[A, -B . Transpose[B]]; // 
      AbsoluteTiming

nrm  =  
        Norm[Xfull . A + Transpose[A] . Xfull + 
            B . Transpose[B]];
If[Head[nrm]  ==  Real,  Print[nrm]]

Solution of Lyapunov equation

Options[LowRankLyapunovSolve]

Xc  =  
          LowRankLyapunovSolve[A, B, SolveMethod -> Krylov,  
            BasisSize -> 50]; // AbsoluteTiming

nrm  =  
        Norm[A . Xc + Xc . Transpose[A] + 
            B . Transpose[B]];
If[Head[nrm]  ==  Real,  Print[nrm]]

Model Reduction

<< LargeScaleControl`RealizationsSparse`
<< NumericalControlMethods`

ss  =  StateSpace[A, B, CC];

Options[DominantSubsystem]

ss1  =  
        DominantSubsystem[ss, 
          ReductionMethod -> LowFactorSquareRoot];

Dimensions[ss1[[1]]]

Bode Plot

tt = BodePlot[Normal[ss]];
tt1 = BodePlot[ss1];

Show[GraphicsArray[{tt, tt1}]]

Riccati Equation

(A^T)X + X A - (XBR^(-1)) (B^T) X  +  (C^T) QC  =  
        0

<< LargeScaleControl`RiccatiSparse`

Q = spdmatgen[1,  {0.1, 20}];
K = N[Table[0. ,  {n*m}, {1}]];
R = N[IdentityMatrix[1]];

X = 
          RiccatiSolve[A,  B, CC,  Q,  R,  K,  
            SolveMethod -> LowRankNewton]; // 
      AbsoluteTiming

Norm[
      Transpose[CC] . Q . CC + Transpose[A] . X + X . A - 
        X . B . R . Transpose[B] . X]